Differential expression of RNA-seq data for a time course analysis on S. mediterranea worms treated with a 70 volt electric field
M. LeGro
4/30/2021
- 1 Filtering and normalizing raw reads
- 2 Plotting sample-to-sample variation
- 3 Limma-voom method for linear modeling of reads
- 4 Row means CPM z-scores for downstream plotting
- 5 Selecting significantly differentially expressed genes using False Discovery Rate (FDR) and Log Fold Change (LogFC) cutoffs
- 6 Creating Heatmaps of all significantly differentially expressed genes for supplement
- 7 Merging pathway annotations with mean z-scores for downstream plotting
- 8 Creating matrices for plotting pathway heatmaps
- 9 Plotting heatmaps of pathways showing significantly differentially expressed transcipts
- 10 Gene set enrichment analysis for 15 minute and 60 minute timepoints
1 Filtering and normalizing raw reads
First I imported the raw counts file and the metadata files and created a DGEList object with the data. Then I filtered the raw counts for all transcripts with a collective CPM value > 0.5 to remove lowly expressing and partial transcripts. The data were normalized to CPM for downstream statistical comparisons.
Attaching package: 'gplots'The following object is masked from 'package:stats':
lowess
Attaching package: 'dplyr'The following objects are masked from 'package:stats':
filter, lagThe following objects are masked from 'package:base':
intersect, setdiff, setequal, unionlibrary(RColorBrewer)
library(dplyr)
setwd("/Volumes/Watsys/EF/")
## Importing raw reads and metadata
counts <- read.delim("2018.8.9_EF_counts.txt", row.names = 1)
meta <- read.delim("2019.11.21_meta.txt", row.names = 1)
meta.f <- factor(meta$Time)
metared <- read.delim("2019.11.16_pDCS_reducedsamps_meta.txt",
row.names = 1)
meta.red <- factor(metared$Time)
## Setting the column names and sample matrix
colnames(counts) <- c("Length", "15min_1", "15min_2", "15min_3",
"30min_1", "30min_2", "30min_3", "60min_1", "60min_2", "60min_3",
"0min_1", "0min_2", "0min_3", "delay_1", "delay_2", "delay_3")
countsmat <- as.matrix(counts[, 2:16])
countsmat15 <- as.matrix(counts[, 2:4])
countsmat60 <- as.matrix(counts[, 8:10])
countsmatdelay <- as.matrix(counts[, 14:16])
countsmatctrl <- as.matrix(counts[, 11:13])
## Creating matrices of samples and a reduced matrix for
## comparisons of 15 and 60 minute timepoints
mat <- cbind(countsmatctrl, countsmat15, countsmat60, countsmatdelay)
colnames(mat) [1] "0min_1" "0min_2" "0min_3" "15min_1" "15min_2" "15min_3" "60min_1"
[8] "60min_2" "60min_3" "delay_1" "delay_2" "delay_3"[1] "0min_1" "0min_2" "0min_3" "15min_1" "15min_2" "15min_3" "60min_1"
[8] "60min_2" "60min_3"# Make an EList object to work with in limma voom
list <- DGEList(mat, group = meta.f) #Creates a DGE list of the counts dataset
listred <- DGEList(matreduced, group = meta.red) #Creates a DGE list of the counts dataset
# Calculate Normalization Factors
list <- calcNormFactors(list)
listred <- calcNormFactors(listred)
# Filter samples with less than 0.5 counts per million in all
# samples
keep <- rowSums(cpm(list) > 0.5) >= 1
list <- list[keep, ]
write.csv(list$counts, file = "2019.11.21_pDCS_filterednormdcounts.csv")
keep <- rowSums(cpm(listred) > 0.5) >= 1
listred <- listred[keep, ]
write.csv(listred$counts, file = "2019.11.16_pDCS_filteredcounts_reducedsamps.csv")
#### Writing the cpm counts into files
CPM <- cpm(list$counts)
write.csv(CPM, file = "2019.11.21_pDCS_countsCPM.csv")
CPMred <- cpm(listred$counts)
write.csv(CPMred, file = "2021.4.15_pDCS_countsCPM_reducedsamps.csv")
countsfiltered <- list$counts
countsfilteredred <- listred$counts
# transpose the matrix so gene rows are columns
CPMt <- t(CPM)
CPMredt <- t(CPMred)
# find the standard score or Z-score of the CPM values for
# downstream graphing and visualization
CPMz <- scale(CPMt, center = TRUE, scale = TRUE)
CPMzscore <- t(CPMz)
CPMredz <- scale(CPMredt, center = TRUE, scale = TRUE)
CPMredzscore <- t(CPMredz)
write.csv(CPMzscore, file = "2021.4.15_pDCS_CPMzscores_allsamps.csv")
write.csv(CPMredzscore, file = "2021.4.15_pDCS_CPMzscores_reducedsamples.csv")2 Plotting sample-to-sample variation
I plotted Multi-dimensional scaling plots of all samples to identify sample-to-sample variation.
cols <- c("black", "red", "green", "blue", "orange")
MDS <- plotMDS(CPM, col = cols[as.numeric(meta.f)])
MDS <- plotMDS(CPMred, col = cols[as.numeric(meta.red)])
MDS <- legend("bottom", legend = levels(meta.f), text.col = cols,
cex = 0.5)
3 Limma-voom method for linear modeling of reads
Next I set up the linear model for statistical comparison based on Chapter 9 of the limma-voom userguide. The Userguides documents can be found here:https://www.bioconductor.org/packages/devel/bioc/vignettes/limma/inst/doc/. I used the ‘voomwithQualityweights’ function which combines the observational-level weights with sample-specific quality weights. This reduces the sample-to-sample variation and downstream removes likelihood of false positives. Bayesian statistics were used to determine the significance of transcripts between samples. After determining the differential expression of transcripts I merged the topTables of statistics with the CPM z-scores.
# Analyzing as for a single factor, section 9.5.2, making
# design matrix
design <- model.matrix(~0 + meta.f)
colnames(design)[1] "meta.fControl" "meta.fdelay" "meta.fmin15" "meta.fmin60" [1] "Control" "delay" "min.15" "min.60" [1] "meta.redControl" "meta.redmin15" "meta.redmin60" [1] "Control" "min.15" "min.60" # Analyzing as for a single factor, section 9.5.2 and
# adjusting for sample-to-sample variation
v <- voomWithQualityWeights(list, design = design, normalize.method = "none",
plot = TRUE)
# fits a linear model to the normalized/filtered dataset
vfit <- lmFit(v, design = design)
vfit <- eBayes(vfit)
vfitred <- lmFit(vred, design = designred)
vfitred <- eBayes(vfitred)
# Setting the statistical contrast matrix
m <- makeContrasts(min15 - Control, min60 - Control, delay -
Control, levels = meta.f)
mred <- makeContrasts(min15 - Control, min60 - Control, levels = meta.red)
# Fitting the statistical contrasts to the linear model
vfit <- contrasts.fit(vfit, contrasts = m)
vfit <- eBayes(vfit)
vfitred <- contrasts.fit(vfitred, contrasts = mred)
vfitred <- eBayes(vfitred)
# test results from testing a set of contrasts equal to 0
results <- decideTests(vfit)
vennDiagram(results)
## printing the t-statistic
tstat <- vfit$t
tstatred <- vfitred$t
# Creating the toptable of test statistics
toptable <- topTable(vfit, number = 3e+05)
toptable_15vctrl <- topTable(vfit, coef = 1, number = 3e+05,
sort.by = "P")
toptable_60vctrl <- topTable(vfit, coef = 2, number = 3e+05,
sort.by = "P")
toptable_delayvctrl <- topTable(vfit, coef = 3, number = 3e+05,
sort.by = "P")
# write.csv(toptable, file = '2018.9.15_EF_toptable.csv')
write.csv(toptable_15vctrl, file = "2019.11.21_pDCS_toptable_coef1_15vctrl.csv")
write.csv(toptable_60vctrl, file = "2019.11.21_pDCS_toptable_coef2_60vctrl.csv")
write.csv(toptable_delayvctrl, file = "2019.11.21_pDCS_toptable_coef3_delayvctrl.csv")
toptablered <- topTable(vfitred, number = 3e+05)
toptable_15vctrlred <- topTable(vfitred, coef = 1, number = 3e+05,
sort.by = "P")
toptable_60vctrlred <- topTable(vfitred, coef = 2, number = 3e+05,
sort.by = "P")
write.csv(toptablered, file = "2021.4.16_pDCS_redsamps_toptable.csv")
write.csv(toptable_15vctrlred, file = "2021.4.16_pDCS_toptable_redsamps_15vctrl.csv")
write.csv(toptable_60vctrlred, file = "2021.4.16_pDCS_toptable_redsamps_60vctrl.csv")
# Merging the topTables of statistics with the CPM z-score
# matrices for all genes
merged15vctrl <- merge.data.frame(toptable_15vctrl, CPMzscore,
by.x = "row.names", by.y = "row.names")
merged60vctrl <- merge.data.frame(toptable_60vctrl, CPMzscore,
by.x = "row.names", by.y = "row.names")
mergeddelayvctrl <- merge.data.frame(toptable_delayvctrl, CPMzscore,
by.x = "row.names", by.y = "row.names")
merged15vctrlred <- merge.data.frame(toptable_15vctrlred, CPMredzscore,
by.x = "row.names", by.y = "row.names")
merged60vctrlred <- merge.data.frame(toptable_60vctrlred, CPMredzscore,
by.x = "row.names", by.y = "row.names")4 Row means CPM z-scores for downstream plotting
The following block of code calculate the mean z-scores for each row of the CPM z-scores matrix. This matrix will be used for downstream plotting of heatmaps.
# calculating the mean z-scores for z-scores counts file
z0 <- rowMeans(subset(CPMredzscore, select = c("0min_1", "0min_2",
"0min_3")), na.rm = TRUE)
z15 <- rowMeans(subset(CPMredzscore, select = c("15min_1", "15min_2",
"15min_3")), na.rm = TRUE)
z60 <- rowMeans(subset(CPMredzscore, select = c("60min_1", "60min_2",
"60min_3")), na.rm = TRUE)
mat <- cbind(z0, z15, z60)
colnames(mat)[1] "z0" "z15" "z60"merged15vctrlred <- merge.data.frame(toptable_15vctrlred, mat,
by.x = "row.names", by.y = "row.names")
merged60vctrlred <- merge.data.frame(toptable_60vctrlred, mat,
by.x = "row.names", by.y = "row.names")
# merged zscores with toptable list
write.csv(merged15vctrlred, file = "2021.4.16_pDCS_meanzscores_15vctrl_CPMmeanzscores_wtoptable.csv")
write.csv(merged60vctrlred, file = "2021.4.16_pDCS_meanzscores_60vctrl_CPMmeanzscores_wtoptable.csv")5 Selecting significantly differentially expressed genes using False Discovery Rate (FDR) and Log Fold Change (LogFC) cutoffs
First the significantly differentially expressed genes are selected using an FDR cutoff <0.05 and <0.01. I decided to use the FDR<0.05 for downstream gene selection since the number of genes is sufficient for downstream plotting. To determine which genes are considered “upregulated” or “downregulated” I’m going to use a logFC cutoff of 0.5 for up or down. However, these cutoffs will only be used for downstream plotting. For the following analysis I used the reduced samples matrices since these time points are needed for publication.
### Finding the numbers of sig DE genes correlating with the
### FDR value of 0.01 and 0.05 for both statistical contrasts
### 1778 sidDE at FDR< 0.05 609 sigDE at FDR< 0.01
ctrlv15_0.05 <- merged15vctrlred[which(merged15vctrlred$adj.P.Val <
0.05), ]
ctrlv15_0.01 <- merged15vctrlred[which(merged15vctrlred$adj.P.Val <
0.01), ]
## 1687 sidDE at FDR< 0.05 sigDE at FDR< 0.01
ctrlv60_0.05 <- merged60vctrlred[which(merged60vctrlred$adj.P.Val <
0.05), ]
ctrlv60_0.01 <- merged60vctrlred[which(merged60vctrlred$adj.P.Val <
0.01), ]
### Finding the upregulated and downregulated genes for each of
### the contrasts and accompanying FDR rates LogFC difference
### of 0.5 and -0.5
# 267 total upregulated genes 315 total downregulated 588
# total updown
ctrlv15_0.05UP_LF0.5 <- ctrlv15_0.05[which(ctrlv15_0.05$logFC >
0.5), ]
ctrlv15_0.05DOWN_LF0.5 <- ctrlv15_0.05[which(ctrlv15_0.05$logFC <
-0.5), ]
ctrlv15_0.05UPDOWN <- subset(ctrlv15_0.05, logFC > 0.5 | logFC <
-0.5)
# 282 total upregulated genes 319 total downregulated 601
# total updown
ctrlv60_0.05UP_LF0.5 <- ctrlv60_0.05[which(ctrlv60_0.05$logFC >
0.5), ]
ctrlv60_0.05DOWN_LF0.5 <- ctrlv60_0.05[which(ctrlv60_0.05$logFC <
-0.5), ]
ctrlv60_0.05UPDOWN <- subset(ctrlv60_0.05, logFC > 0.5 | logFC <
-0.5)6 Creating Heatmaps of all significantly differentially expressed genes for supplement
Next I created matrices for plotting heatmaps of all significantly differentially expressed genes with an FDR <0.05. I also created complimentary heatmaps with a LogFC cutoff of 0.5 up and down.
### Creating the matrices
### ##########################################################
ctrlv15_0.05allhm <- cbind(ctrlv15_0.05$z0, ctrlv15_0.05$z15)
rownames(ctrlv15_0.05allhm) <- c(ctrlv15_0.05$Row.names)
colnames(ctrlv15_0.05allhm) <- c("Control", "15_minute")
ctrlv15_0.05UPDOWNhm <- cbind(ctrlv15_0.05UPDOWN$z0, ctrlv15_0.05UPDOWN$z15)
rownames(ctrlv15_0.05UPDOWNhm) <- c(ctrlv15_0.05UPDOWN$Row.names)
colnames(ctrlv15_0.05UPDOWNhm) <- c("Control", "15_minute")
ctrlv60_0.05allhm <- cbind(ctrlv60_0.05$z0, ctrlv60_0.05$z60)
rownames(ctrlv60_0.05allhm) <- c(ctrlv60_0.05$Row.names)
colnames(ctrlv60_0.05allhm) <- c("Control", "60_minute")
ctrlv60_0.05UPDOWNhm <- cbind(ctrlv60_0.05UPDOWN$z0, ctrlv60_0.05UPDOWN$z60)
rownames(ctrlv60_0.05UPDOWNhm) <- c(ctrlv60_0.05UPDOWN$Row.names)
colnames(ctrlv60_0.05UPDOWNhm) <- c("Control", "60_minute")
### Creating the heatmaps
### ##########################################################
### Calling the packages and setting the colors
library(ComplexHeatmap)Loading required package: grid========================================
ComplexHeatmap version 2.6.2
Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
Github page: https://github.com/jokergoo/ComplexHeatmap
Documentation: http://jokergoo.github.io/ComplexHeatmap-reference
If you use it in published research, please cite:
Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional
genomic data. Bioinformatics 2016.
This message can be suppressed by:
suppressPackageStartupMessages(library(ComplexHeatmap))
================================================================================
circlize version 0.4.12
CRAN page: https://cran.r-project.org/package=circlize
Github page: https://github.com/jokergoo/circlize
Documentation: https://jokergoo.github.io/circlize_book/book/
If you use it in published research, please cite:
Gu, Z. circlize implements and enhances circular visualization
in R. Bioinformatics 2014.
This message can be suppressed by:
suppressPackageStartupMessages(library(circlize))
========================================library(RColorBrewer)
mycolz <- colorRamp2(c(-1, 0, 1), c("purple", "white", "dark cyan"))
##### Heatmaps of sigDE genes for each comparison plotting the
##### mean zscores####
Heatmap(ctrlv15_0.05allhm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(ctrlv15_0.05UPDOWNhm, name = "Z-score", col = mycolz,
column_names_gp = gpar(fontsize = 8), cluster_columns = FALSE,
cluster_rows = TRUE)
Heatmap(ctrlv60_0.05allhm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(ctrlv60_0.05UPDOWNhm, name = "Z-score", col = mycolz,
column_names_gp = gpar(fontsize = 8), cluster_columns = FALSE,
cluster_rows = TRUE)
7 Merging pathway annotations with mean z-scores for downstream plotting
The following lists of pathway annotations are Homo sapien ortholog BLAST annotations that were downloaded from the planmine database here:https://www.bioconductor.org/packages/devel/bioc/vignettes/limma/inst/doc/. I merged the pathway annotations with the significantly differentially expressed transcripts for the 15 minute and 60 minute pairwise comparisons (FDR<0.05 and additional lists with LogFC cutoff of 0.5).
celldeath <- read.delim("2020.8.26_celldeath_ddsmedv6_homosapiens_blastannots.txt",
row.names = 1)
neural <- read.delim("2020.8.26_neural_ddsmedv6_homosapiens_blastannots.txt",
row.names = 1)
proliferation <- read.delim("2020.8.26_proliferation_ddsmedv6_homosapiens_blastannots.txt",
row.names = 1)
replication <- read.delim("2020.8.26_replication_ddsmedv6_homosapiens_blastannots.txt",
row.names = 1)
signaling <- read.delim("2020.8.26_signaling_ddsmedv6_homosapiens_blastannots.txt",
row.names = 1)
calcium <- read.delim("2020.11.2_calcium_v6homologs_homosapiens_blastannots.txt",
row.names = 1)
cellmigration <- read.delim("2020.11.2_cellmigration_v6homologs_homosapiens_blastannots.txt",
row.names = 1)
DNAdamage <- read.delim("2020.11.2_DNAdamage_v6homologs_homosapiens_blastannots.txt",
row.names = 1)
nbsc <- read.delim("2020.11.2_tspanpaper_markersfulllist_NBSConly_v6homologs.txt",
row.names = 1)
sublethal <- read.delim("2020.11.2_tspanpaper_markersfulllist_SLonly_v6homologs.txt",
row.names = 1)
cellcycle <- read.delim("2020.11.16_cellcycle_geneIDs_homosapiens_blastannots.txt",
row.names = 1)
dedifferentiation <- read.delim("2020.11.16_de-differentiation_geneIDs_homosapiens_blastannots.txt",
row.names = 1)
DNArepair <- read.delim("2020.11.16_DNArepair_geneIDs_homosapiens_blastannots.txt",
row.names = 1)
HR <- read.delim("2020.11.16_homologousrecombination_geneIDs_homosapiens_blastannots.txt",
row.names = 1)
MMR <- read.delim("2020.11.16_mismatchrepair_geneIDs_homosapiens_blastannots.txt",
row.names = 1)
ieg <- read.delim("2020.11.16_pDCS_ieglist.txt", row.names = 1)
piwi <- read.delim("2020.11.16_piwi_geneIDs_blastannots.txt",
row.names = 1)
tspan <- read.delim("2020.11.16_t-span_geneIDs_blastannots.txt",
row.names = 1)
tetraspanin <- read.delim("2020.11.16_tetraspanin_geneIDs_homosapiens_blastannots.txt",
row.names = 1)
NHEJ <- read.delim("2020.11.16_NHEJ_geneIDs_homosapiens_blastannots.txt",
row.names = 1)
mex3 <- read.delim("2020.11.22_EF_Hippo_mex3b.txt", row.names = 1)
cellmigS6 <- read.delim("2020.11.29_pDCS_cellmigration_figureS6_homologs_blastannots.txt",
row.names = 1)
rad51 <- read.delim("2020.11.29_pDCS_rad51_homologs_blastannots.txt",
row.names = 1)
agat <- read.delim("2021.2.15_agat_blastannots.txt", row.names = 1)
NBFIG <- read.delim("2021.2.22_stemcellmarkers.txt", row.names = 1)
#### MERGING THE gene homologs for pathways lists with the sigDE
#### genes for each timepoint and the specifically upregulated
#### and downregulated ##### cell death ####
celldeath15 <- merge.data.frame(ctrlv15_0.05, celldeath, by.x = "Row.names",
by.y = "row.names")
celldeath60 <- merge.data.frame(ctrlv60_0.05, celldeath, by.x = "Row.names",
by.y = "row.names")
celldeath15UPDOWN <- merge.data.frame(ctrlv15_0.05UPDOWN, celldeath,
by.x = "Row.names", by.y = "row.names")
celldeath60UPDOWN <- merge.data.frame(ctrlv60_0.05UPDOWN, celldeath,
by.x = "Row.names", by.y = "row.names")
#### neural ####
neural15 <- merge.data.frame(ctrlv15_0.05, neural, by.x = "Row.names",
by.y = "row.names")
neural60 <- merge.data.frame(ctrlv60_0.05, neural, by.x = "Row.names",
by.y = "row.names")
neural15UPDOWN <- merge.data.frame(ctrlv15_0.05UPDOWN, neural,
by.x = "Row.names", by.y = "row.names")
neural60UPDOWN <- merge.data.frame(ctrlv60_0.05UPDOWN, neural,
by.x = "Row.names", by.y = "row.names")
#### proliferation ####
proliferation15 <- merge.data.frame(ctrlv15_0.05, proliferation,
by.x = "Row.names", by.y = "row.names")
proliferation60 <- merge.data.frame(ctrlv60_0.05, proliferation,
by.x = "Row.names", by.y = "row.names")
proliferation15UPDOWN <- merge.data.frame(ctrlv15_0.05UPDOWN,
proliferation, by.x = "Row.names", by.y = "row.names")
proliferation60UPDOWN <- merge.data.frame(ctrlv60_0.05UPDOWN,
proliferation, by.x = "Row.names", by.y = "row.names")
#### replication ####
replication15 <- merge.data.frame(ctrlv15_0.05, replication,
by.x = "Row.names", by.y = "row.names")
replication60 <- merge.data.frame(ctrlv60_0.05, replication,
by.x = "Row.names", by.y = "row.names")
replication15UPDOWN <- merge.data.frame(ctrlv15_0.05UPDOWN, replication,
by.x = "Row.names", by.y = "row.names")
replication60UPDOWN <- merge.data.frame(ctrlv60_0.05UPDOWN, replication,
by.x = "Row.names", by.y = "row.names")
#### signaling ####
signaling15 <- merge.data.frame(ctrlv15_0.05, signaling, by.x = "Row.names",
by.y = "row.names")
signaling60 <- merge.data.frame(ctrlv60_0.05, signaling, by.x = "Row.names",
by.y = "row.names")
signaling15UPDOWN <- merge.data.frame(ctrlv15_0.05UPDOWN, signaling,
by.x = "Row.names", by.y = "row.names")
signaling60UPDOWN <- merge.data.frame(ctrlv60_0.05UPDOWN, signaling,
by.x = "Row.names", by.y = "row.names")
#### calcium ####
calcium15 <- merge.data.frame(ctrlv15_0.05, calcium, by.x = "Row.names",
by.y = "row.names")
calcium60 <- merge.data.frame(ctrlv60_0.05, calcium, by.x = "Row.names",
by.y = "row.names")
calcium15UPDOWN <- merge.data.frame(ctrlv15_0.05UPDOWN, calcium,
by.x = "Row.names", by.y = "row.names")
calcium60UPDOWN <- merge.data.frame(ctrlv60_0.05UPDOWN, calcium,
by.x = "Row.names", by.y = "row.names")
#### cell migration ####
cellmigration15 <- merge.data.frame(ctrlv15_0.05, cellmigration,
by.x = "Row.names", by.y = "row.names")
cellmigration60 <- merge.data.frame(ctrlv60_0.05, cellmigration,
by.x = "Row.names", by.y = "row.names")
cellmigration15UPDOWN <- merge.data.frame(ctrlv15_0.05UPDOWN,
cellmigration, by.x = "Row.names", by.y = "row.names")
cellmigration60UPDOWN <- merge.data.frame(ctrlv60_0.05UPDOWN,
cellmigration, by.x = "Row.names", by.y = "row.names")
#### DNA damage ####
DNAdamage15 <- merge.data.frame(ctrlv15_0.05, DNAdamage, by.x = "Row.names",
by.y = "row.names")
DNAdamage60 <- merge.data.frame(ctrlv60_0.05, DNAdamage, by.x = "Row.names",
by.y = "row.names")
DNAdamage15UPDOWN <- merge.data.frame(ctrlv15_0.05UPDOWN, DNAdamage,
by.x = "Row.names", by.y = "row.names")
DNAdamage60UPDOWN <- merge.data.frame(ctrlv60_0.05UPDOWN, DNAdamage,
by.x = "Row.names", by.y = "row.names")
#### cellcycle ####
cellcycle15 <- merge.data.frame(ctrlv15_0.05, cellcycle, by.x = "Row.names",
by.y = "row.names")
cellcycle60 <- merge.data.frame(ctrlv60_0.05, cellcycle, by.x = "Row.names",
by.y = "row.names")
cellcycle15UPDOWN <- merge.data.frame(ctrlv15_0.05UPDOWN, cellcycle,
by.x = "Row.names", by.y = "row.names")
cellcycle60UPDOWN <- merge.data.frame(ctrlv60_0.05UPDOWN, cellcycle,
by.x = "Row.names", by.y = "row.names")
#### DNArepair ####
DNArepair15 <- merge.data.frame(ctrlv15_0.05, DNArepair, by.x = "Row.names",
by.y = "row.names")
DNArepair60 <- merge.data.frame(ctrlv60_0.05, DNArepair, by.x = "Row.names",
by.y = "row.names")
DNArepair15UPDOWN <- merge.data.frame(ctrlv15_0.05UPDOWN, DNArepair,
by.x = "Row.names", by.y = "row.names")
DNArepair60UPDOWN <- merge.data.frame(ctrlv60_0.05UPDOWN, DNArepair,
by.x = "Row.names", by.y = "row.names")
#### HR ####
HR15 <- merge.data.frame(ctrlv15_0.05, HR, by.x = "Row.names",
by.y = "row.names")
HR60 <- merge.data.frame(ctrlv60_0.05, HR, by.x = "Row.names",
by.y = "row.names")
HR15UPDOWN <- merge.data.frame(ctrlv15_0.05UPDOWN, HR, by.x = "Row.names",
by.y = "row.names")
HR60UPDOWN <- merge.data.frame(ctrlv60_0.05UPDOWN, HR, by.x = "Row.names",
by.y = "row.names")
#### MMR ####
MMR15 <- merge.data.frame(ctrlv15_0.05, MMR, by.x = "Row.names",
by.y = "row.names")
MMR60 <- merge.data.frame(ctrlv60_0.05, MMR, by.x = "Row.names",
by.y = "row.names")
MMR15UPDOWN <- merge.data.frame(ctrlv15_0.05UPDOWN, MMR, by.x = "Row.names",
by.y = "row.names")
MMR60UPDOWN <- merge.data.frame(ctrlv60_0.05UPDOWN, MMR, by.x = "Row.names",
by.y = "row.names")
#### NHEJ ####
NHEJ15 <- merge.data.frame(ctrlv15_0.05, NHEJ, by.x = "Row.names",
by.y = "row.names")
NHEJ60 <- merge.data.frame(ctrlv60_0.05, NHEJ, by.x = "Row.names",
by.y = "row.names")
NHEJ15UPDOWN <- merge.data.frame(ctrlv15_0.05UPDOWN, NHEJ, by.x = "Row.names",
by.y = "row.names")
NHEJ60UPDOWN <- merge.data.frame(ctrlv60_0.05UPDOWN, NHEJ, by.x = "Row.names",
by.y = "row.names")
#### cell migration figure S6 ####
figurecellmig15 <- merge.data.frame(ctrlv15_0.05, cellmigS6,
by.x = "Row.names", by.y = "row.names")
figurecellmig60 <- merge.data.frame(ctrlv60_0.05, cellmigS6,
by.x = "Row.names", by.y = "row.names")
figurecellmig15UPDOWN <- merge.data.frame(ctrlv15_0.05UPDOWN,
cellmigS6, by.x = "Row.names", by.y = "row.names")
figurecellmig60UPDOWN <- merge.data.frame(ctrlv60_0.05UPDOWN,
cellmigS6, by.x = "Row.names", by.y = "row.names")
#### NB markers figure 1 ####
NBFIG15 <- merge.data.frame(ctrlv15_0.05, NBFIG, by.x = "Row.names",
by.y = "row.names")
NBFIG60 <- merge.data.frame(ctrlv60_0.05, NBFIG, by.x = "Row.names",
by.y = "row.names")
NBFIG15UPDOWN <- merge.data.frame(ctrlv15_0.05UPDOWN, NBFIG,
by.x = "Row.names", by.y = "row.names")
NBFIG60UPDOWN <- merge.data.frame(ctrlv60_0.05UPDOWN, NBFIG,
by.x = "Row.names", by.y = "row.names")
#### rad51 ####
rad5115 <- merge.data.frame(ctrlv15_0.05, rad51, by.x = "Row.names",
by.y = "row.names")
rad5160 <- merge.data.frame(ctrlv60_0.05, rad51, by.x = "Row.names",
by.y = "row.names")
rad5115UPDOWN <- merge.data.frame(ctrlv15_0.05UPDOWN, rad51,
by.x = "Row.names", by.y = "row.names")
rad5160UPDOWN <- merge.data.frame(ctrlv60_0.05UPDOWN, rad51,
by.x = "Row.names", by.y = "row.names")
#### nbsc ####
nbsc15 <- merge.data.frame(ctrlv15_0.05, nbsc, by.x = "Row.names",
by.y = "row.names")
nbsc60 <- merge.data.frame(ctrlv60_0.05, nbsc, by.x = "Row.names",
by.y = "row.names")
nbsc15UPDOWN <- merge.data.frame(ctrlv15_0.05UPDOWN, nbsc, by.x = "Row.names",
by.y = "row.names")
nbsc60UPDOWN <- merge.data.frame(ctrlv60_0.05UPDOWN, nbsc, by.x = "Row.names",
by.y = "row.names")
#### sublethal ####
sublethal15 <- merge.data.frame(ctrlv15_0.05, sublethal, by.x = "Row.names",
by.y = "row.names")
sublethal60 <- merge.data.frame(ctrlv60_0.05, sublethal, by.x = "Row.names",
by.y = "row.names")
sublethal15UPDOWN <- merge.data.frame(ctrlv15_0.05UPDOWN, sublethal,
by.x = "Row.names", by.y = "row.names")
sublethal60UPDOWN <- merge.data.frame(ctrlv60_0.05UPDOWN, sublethal,
by.x = "Row.names", by.y = "row.names")
#### ieg ####
ieg15 <- merge.data.frame(ctrlv15_0.05, ieg, by.x = "Row.names",
by.y = "row.names")
ieg60 <- merge.data.frame(ctrlv60_0.05, ieg, by.x = "Row.names",
by.y = "row.names")
ieg15UPDOWN <- merge.data.frame(ctrlv15_0.05UPDOWN, ieg, by.x = "Row.names",
by.y = "row.names")
ieg60UPDOWN <- merge.data.frame(ctrlv60_0.05UPDOWN, ieg, by.x = "Row.names",
by.y = "row.names")
#### piwi ####
piwi15 <- merge.data.frame(ctrlv15_0.05, piwi, by.x = "Row.names",
by.y = "row.names")
piwi60 <- merge.data.frame(ctrlv60_0.05, piwi, by.x = "Row.names",
by.y = "row.names")
piwi15UPDOWN <- merge.data.frame(ctrlv15_0.05UPDOWN, piwi, by.x = "Row.names",
by.y = "row.names")
piwi60UPDOWN <- merge.data.frame(ctrlv60_0.05UPDOWN, piwi, by.x = "Row.names",
by.y = "row.names")8 Creating matrices for plotting pathway heatmaps
The following chunk of code creates matrices for 15 and 60 minute contrasts using all significantly differentially expressed genes and additional heatmaps using a LogFC cutoff of 0.5.
### Creating the matrices
### ##########################################################
### Cell death
### ##################################################
celldeath15$descrip <- paste(celldeath15$blast_description, celldeath15$Row.names,
sep = "_")
celldeath15$symbol <- paste(celldeath15$blast_symbol, celldeath15$Row.names,
sep = "_")
celldeath15_hm <- cbind(celldeath15$z0, celldeath15$z15)
rownames(celldeath15_hm) <- c(celldeath15$symbol)
colnames(celldeath15_hm) <- c("Control", "15_minute")
celldeath60$descrip <- paste(celldeath60$blast_description, celldeath60$Row.names,
sep = "_")
celldeath60$symbol <- paste(celldeath60$blast_symbol, celldeath60$Row.names,
sep = "_")
celldeath60_hm <- cbind(celldeath60$z0, celldeath60$z15)
rownames(celldeath60_hm) <- c(celldeath60$symbol)
colnames(celldeath60_hm) <- c("Control", "60_minute")
celldeath15UPDOWN$descrip <- paste(celldeath15UPDOWN$blast_description,
celldeath15UPDOWN$Row.names, sep = "_")
celldeath15UPDOWN$symbol <- paste(celldeath15UPDOWN$blast_symbol,
celldeath15UPDOWN$Row.names, sep = "_")
celldeath15UPDOWN_hm <- cbind(celldeath15UPDOWN$z0, celldeath15UPDOWN$z15)
rownames(celldeath15UPDOWN_hm) <- c(celldeath15UPDOWN$symbol)
colnames(celldeath15UPDOWN_hm) <- c("Control", "15_minute")
celldeath60UPDOWN$descrip <- paste(celldeath60UPDOWN$blast_description,
celldeath60UPDOWN$Row.names, sep = "_")
celldeath60UPDOWN$symbol <- paste(celldeath60UPDOWN$blast_symbol,
celldeath60UPDOWN$Row.names, sep = "_")
celldeath60UPDOWN_hm <- cbind(celldeath60UPDOWN$z0, celldeath60UPDOWN$z15)
rownames(celldeath60UPDOWN_hm) <- c(celldeath60UPDOWN$symbol)
colnames(celldeath60UPDOWN_hm) <- c("Control", "60_minute")
################################################################################ neural ##################################################
neural15_hm <- cbind(neural15$z0, neural15$z15)
rownames(neural15_hm) <- c(neural15$Blast_Symbol)
colnames(neural15_hm) <- c("Control", "15_minute")
neural60_hm <- cbind(neural60$z0, neural60$z15)
rownames(neural60_hm) <- c(neural60$Blast_Symbol)
colnames(neural60_hm) <- c("Control", "60_minute")
neural15UPDOWN_hm <- cbind(neural15UPDOWN$z0, neural15UPDOWN$z15)
rownames(neural15UPDOWN_hm) <- c(neural15UPDOWN$Blast_Symbol)
colnames(neural15UPDOWN_hm) <- c("Control", "15_minute")
neural60UPDOWN_hm <- cbind(neural60UPDOWN$z0, neural60UPDOWN$z15)
rownames(neural60UPDOWN_hm) <- c(neural60UPDOWN$Blast_Symbol)
colnames(neural60UPDOWN_hm) <- c("Control", "60_minute")
################################################################################ proliferation
################################################################################ ##################################################3
proliferation15$descrip <- paste(proliferation15$blast_description,
proliferation15$Row.names, sep = "_")
proliferation15$symbol <- paste(proliferation15$blast_symbol,
proliferation15$Row.names, sep = "_")
proliferation15_hm <- cbind(proliferation15$z0, proliferation15$z15)
rownames(proliferation15_hm) <- c(proliferation15$symbol)
colnames(proliferation15_hm) <- c("Control", "15_minute")
proliferation60$descrip <- paste(proliferation60$blast_description,
proliferation60$Row.names, sep = "_")
proliferation60$symbol <- paste(proliferation60$blast_symbol,
proliferation60$Row.names, sep = "_")
proliferation60_hm <- cbind(proliferation60$z0, proliferation60$z15)
rownames(proliferation60_hm) <- c(proliferation60$symbol)
colnames(proliferation60_hm) <- c("Control", "60_minute")
proliferation15UPDOWN$descrip <- paste(proliferation15UPDOWN$blast_description,
proliferation15UPDOWN$Row.names, sep = "_")
proliferation15UPDOWN$symbol <- paste(proliferation15UPDOWN$blast_symbol,
proliferation15UPDOWN$Row.names, sep = "_")
proliferation15UPDOWN_hm <- cbind(proliferation15UPDOWN$z0, proliferation15UPDOWN$z15)
rownames(proliferation15UPDOWN_hm) <- c(proliferation15UPDOWN$symbol)
colnames(proliferation15UPDOWN_hm) <- c("Control", "15_minute")
proliferation60UPDOWN$descrip <- paste(proliferation60UPDOWN$blast_description,
proliferation60UPDOWN$Row.names, sep = "_")
proliferation60UPDOWN$symbol <- paste(proliferation60UPDOWN$blast_symbol,
proliferation60UPDOWN$Row.names, sep = "_")
proliferation60UPDOWN_hm <- cbind(proliferation60UPDOWN$z0, proliferation60UPDOWN$z15)
rownames(proliferation60UPDOWN_hm) <- c(proliferation60UPDOWN$symbol)
colnames(proliferation60UPDOWN_hm) <- c("Control", "60_minute")
################################################################################ replication
################################################################################ ##################################################3
replication15_hm <- cbind(replication15$z0, replication15$z15)
rownames(replication15_hm) <- c(replication15$Blast_Symbol)
colnames(replication15_hm) <- c("Control", "15_minute")
replication60_hm <- cbind(replication60$z0, replication60$z15)
rownames(replication60_hm) <- c(replication60$Blast_Symbol)
colnames(replication60_hm) <- c("Control", "60_minute")
replication15UPDOWN_hm <- cbind(replication15UPDOWN$z0, replication15UPDOWN$z15)
rownames(replication15UPDOWN_hm) <- c(replication15UPDOWN$Blast_Symbol)
colnames(replication15UPDOWN_hm) <- c("Control", "15_minute")
replication60UPDOWN_hm <- cbind(replication60UPDOWN$z0, replication60UPDOWN$z15)
rownames(replication60UPDOWN_hm) <- c(replication60UPDOWN$Blast_Symbol)
colnames(replication60UPDOWN_hm) <- c("Control", "60_minute")
################################################################################ signaling
################################################################################ ##################################################3
signaling15_hm <- cbind(signaling15$z0, signaling15$z15)
rownames(signaling15_hm) <- c(signaling15$Blast_Symbol)
colnames(signaling15_hm) <- c("Control", "15_minute")
signaling60_hm <- cbind(signaling60$z0, signaling60$z15)
rownames(signaling60_hm) <- c(signaling60$Blast_Symbol)
colnames(signaling60_hm) <- c("Control", "60_minute")
signaling15UPDOWN_hm <- cbind(signaling15UPDOWN$z0, signaling15UPDOWN$z15)
rownames(signaling15UPDOWN_hm) <- c(signaling15UPDOWN$Blast_Symbol)
colnames(signaling15UPDOWN_hm) <- c("Control", "15_minute")
signaling60UPDOWN_hm <- cbind(signaling60UPDOWN$z0, signaling60UPDOWN$z15)
rownames(signaling60UPDOWN_hm) <- c(signaling60UPDOWN$Blast_Symbol)
colnames(signaling60UPDOWN_hm) <- c("Control", "60_minute")
################################################################################ calcium ##################################################3
calcium15_hm <- cbind(calcium15$z0, calcium15$z15)
rownames(calcium15_hm) <- c(calcium15$Blast_Symbol)
colnames(calcium15_hm) <- c("Control", "15_minute")
calcium60_hm <- cbind(calcium60$z0, calcium60$z15)
rownames(calcium60_hm) <- c(calcium60$Blast_Symbol)
colnames(calcium60_hm) <- c("Control", "60_minute")
calcium15UPDOWN_hm <- cbind(calcium15UPDOWN$z0, calcium15UPDOWN$z15)
rownames(calcium15UPDOWN_hm) <- c(calcium15UPDOWN$Blast_Symbol)
colnames(calcium15UPDOWN_hm) <- c("Control", "15_minute")
calcium60UPDOWN_hm <- cbind(calcium60UPDOWN$z0, calcium60UPDOWN$z15)
rownames(calcium60UPDOWN_hm) <- c(calcium60UPDOWN$Blast_Symbol)
colnames(calcium60UPDOWN_hm) <- c("Control", "60_minute")
################################################################################ cellmigration
################################################################################ ##################################################3
cellmigration15_hm <- cbind(cellmigration15$z0, cellmigration15$z15)
rownames(cellmigration15_hm) <- c(cellmigration15$Blast_Symbol)
colnames(cellmigration15_hm) <- c("Control", "15_minute")
cellmigration60_hm <- cbind(cellmigration60$z0, cellmigration60$z15)
rownames(cellmigration60_hm) <- c(cellmigration60$Blast_Symbol)
colnames(cellmigration60_hm) <- c("Control", "60_minute")
cellmigration15UPDOWN_hm <- cbind(cellmigration15UPDOWN$z0, cellmigration15UPDOWN$z15)
rownames(cellmigration15UPDOWN_hm) <- c(cellmigration15UPDOWN$Blast_Symbol)
colnames(cellmigration15UPDOWN_hm) <- c("Control", "15_minute")
cellmigration60UPDOWN_hm <- cbind(cellmigration60UPDOWN$z0, cellmigration60UPDOWN$z15)
rownames(cellmigration60UPDOWN_hm) <- c(cellmigration60UPDOWN$Blast_Symbol)
colnames(cellmigration60UPDOWN_hm) <- c("Control", "60_minute")
################################################################################ DNAdamage
################################################################################ ##################################################3
DNAdamage15_hm <- cbind(DNAdamage15$z0, DNAdamage15$z15)
rownames(DNAdamage15_hm) <- c(DNAdamage15$Blast_Symbol)
colnames(DNAdamage15_hm) <- c("Control", "15_minute")
DNAdamage60_hm <- cbind(DNAdamage60$z0, DNAdamage60$z15)
rownames(DNAdamage60_hm) <- c(DNAdamage60$Blast_Symbol)
colnames(DNAdamage60_hm) <- c("Control", "60_minute")
DNAdamage15UPDOWN_hm <- cbind(DNAdamage15UPDOWN$z0, DNAdamage15UPDOWN$z15)
rownames(DNAdamage15UPDOWN_hm) <- c(DNAdamage15UPDOWN$Blast_Symbol)
colnames(DNAdamage15UPDOWN_hm) <- c("Control", "15_minute")
DNAdamage60UPDOWN_hm <- cbind(DNAdamage60UPDOWN$z0, DNAdamage60UPDOWN$z15)
rownames(DNAdamage60UPDOWN_hm) <- c(DNAdamage60UPDOWN$Blast_Symbol)
colnames(DNAdamage60UPDOWN_hm) <- c("Control", "60_minute")
################################################################################ nbsc ##################################################3
nbsc15_hm <- cbind(nbsc15$z0, nbsc15$z15)
rownames(nbsc15_hm) <- c(nbsc15$Row.names)
colnames(nbsc15_hm) <- c("Control", "15_minute")
nbsc60_hm <- cbind(nbsc60$z0, nbsc60$z15)
rownames(nbsc60_hm) <- c(nbsc60$Row.names)
colnames(nbsc60_hm) <- c("Control", "60_minute")
nbsc15UPDOWN_hm <- cbind(nbsc15UPDOWN$z0, nbsc15UPDOWN$z15)
rownames(nbsc15UPDOWN_hm) <- c(nbsc15UPDOWN$Row.names)
colnames(nbsc15UPDOWN_hm) <- c("Control", "15_minute")
nbsc60UPDOWN_hm <- cbind(nbsc60UPDOWN$z0, nbsc60UPDOWN$z15)
rownames(nbsc60UPDOWN_hm) <- c(nbsc60UPDOWN$Row.names)
colnames(nbsc60UPDOWN_hm) <- c("Control", "60_minute")
################################################################################ sublethal
################################################################################ ##################################################3
sublethal15_hm <- cbind(sublethal15$z0, sublethal15$z15)
rownames(sublethal15_hm) <- c(sublethal15$Row.names)
colnames(sublethal15_hm) <- c("Control", "15_minute")
sublethal60_hm <- cbind(sublethal60$z0, sublethal60$z15)
rownames(sublethal60_hm) <- c(sublethal60$Row.names)
colnames(sublethal60_hm) <- c("Control", "60_minute")
sublethal15UPDOWN_hm <- cbind(sublethal15UPDOWN$z0, sublethal15UPDOWN$z15)
rownames(sublethal15UPDOWN_hm) <- c(sublethal15UPDOWN$Row.names)
colnames(sublethal15UPDOWN_hm) <- c("Control", "15_minute")
sublethal60UPDOWN_hm <- cbind(sublethal60UPDOWN$z0, sublethal60UPDOWN$z15)
rownames(sublethal60UPDOWN_hm) <- c(sublethal60UPDOWN$Row.names)
colnames(sublethal60UPDOWN_hm) <- c("Control", "60_minute")
################################################################################ cellcycle
################################################################################ ##################################################3
cellcycle15_hm <- cbind(cellcycle15$z0, cellcycle15$z15)
rownames(cellcycle15_hm) <- c(cellcycle15$Transcript...Blast.Sequence.Features...Blast.Domain...Symbol)
colnames(cellcycle15_hm) <- c("Control", "15_minute")
cellcycle60_hm <- cbind(cellcycle60$z0, cellcycle60$z15)
rownames(cellcycle60_hm) <- c(cellcycle60$Transcript...Blast.Sequence.Features...Blast.Domain...Symbol)
colnames(cellcycle60_hm) <- c("Control", "60_minute")
cellcycle15UPDOWN_hm <- cbind(cellcycle15UPDOWN$z0, cellcycle15UPDOWN$z15)
rownames(cellcycle15UPDOWN_hm) <- c(cellcycle15UPDOWN$Transcript...Blast.Sequence.Features...Blast.Domain...Symbol)
colnames(cellcycle15UPDOWN_hm) <- c("Control", "15_minute")
cellcycle60UPDOWN_hm <- cbind(cellcycle60UPDOWN$z0, cellcycle60UPDOWN$z15)
rownames(cellcycle60UPDOWN_hm) <- c(cellcycle60UPDOWN$Transcript...Blast.Sequence.Features...Blast.Domain...Symbol)
colnames(cellcycle60UPDOWN_hm) <- c("Control", "60_minute")
################################################################################ DNArepair
################################################################################ ##################################################3
DNArepair15_hm <- cbind(DNArepair15$z0, DNArepair15$z15)
rownames(DNArepair15_hm) <- c(DNArepair15$Transcript...Blast.Sequence.Features...Blast.Domain...Symbol)
colnames(DNArepair15_hm) <- c("Control", "15_minute")
DNArepair60_hm <- cbind(DNArepair60$z0, DNArepair60$z15)
rownames(DNArepair60_hm) <- c(DNArepair60$Transcript...Blast.Sequence.Features...Blast.Domain...Symbol)
colnames(DNArepair60_hm) <- c("Control", "60_minute")
DNArepair15UPDOWN_hm <- cbind(DNArepair15UPDOWN$z0, DNArepair15UPDOWN$z15)
rownames(DNArepair15UPDOWN_hm) <- c(DNArepair15UPDOWN$Transcript...Blast.Sequence.Features...Blast.Domain...Symbol)
colnames(DNArepair15UPDOWN_hm) <- c("Control", "15_minute")
DNArepair60UPDOWN_hm <- cbind(DNArepair60UPDOWN$z0, DNArepair60UPDOWN$z15)
rownames(DNArepair60UPDOWN_hm) <- c(DNArepair60UPDOWN$Transcript...Blast.Sequence.Features...Blast.Domain...Symbol)
colnames(DNArepair60UPDOWN_hm) <- c("Control", "60_minute")
################################################################################ HR ##################################################3
HR15_hm <- cbind(HR15$z0, HR15$z15)
rownames(HR15_hm) <- c(HR15$Transcript...Blast.Sequence.Features...Blast.Domain...Symbol)
colnames(HR15_hm) <- c("Control", "15_minute")
HR60_hm <- cbind(HR60$z0, HR60$z15)
rownames(HR60_hm) <- c(HR60$Transcript...Blast.Sequence.Features...Blast.Domain...Symbol)
colnames(HR60_hm) <- c("Control", "60_minute")
HR15UPDOWN_hm <- cbind(HR15UPDOWN$z0, HR15UPDOWN$z15)
rownames(HR15UPDOWN_hm) <- c(HR15UPDOWN$Transcript...Blast.Sequence.Features...Blast.Domain...Symbol)
colnames(HR15UPDOWN_hm) <- c("Control", "15_minute")
HR60UPDOWN_hm <- cbind(HR60UPDOWN$z0, HR60UPDOWN$z15)
rownames(HR60UPDOWN_hm) <- c(HR60UPDOWN$Transcript...Blast.Sequence.Features...Blast.Domain...Symbol)
colnames(HR60UPDOWN_hm) <- c("Control", "60_minute")
################################################################################ MMR ##################################################3
MMR15_hm <- cbind(MMR15$z0, MMR15$z15)
rownames(MMR15_hm) <- c(MMR15$Transcript...Blast.Sequence.Features...Blast.Domain...Symbol)
colnames(MMR15_hm) <- c("Control", "15_minute")
MMR60_hm <- cbind(MMR60$z0, MMR60$z15)
rownames(MMR60_hm) <- c(MMR60$Transcript...Blast.Sequence.Features...Blast.Domain...Symbol)
colnames(MMR60_hm) <- c("Control", "60_minute")
MMR15UPDOWN_hm <- cbind(MMR15UPDOWN$z0, MMR15UPDOWN$z15)
rownames(MMR15UPDOWN_hm) <- c(MMR15UPDOWN$Transcript...Blast.Sequence.Features...Blast.Domain...Symbol)
colnames(MMR15UPDOWN_hm) <- c("Control", "15_minute")
MMR60UPDOWN_hm <- cbind(MMR60UPDOWN$z0, MMR60UPDOWN$z15)
rownames(MMR60UPDOWN_hm) <- c(MMR60UPDOWN$Transcript...Blast.Sequence.Features...Blast.Domain...Symbol)
colnames(MMR60UPDOWN_hm) <- c("Control", "60_minute")
################################################################################ ieg ##################################################3
ieg15_hm <- cbind(ieg15$z0, ieg15$z15)
rownames(ieg15_hm) <- c(ieg15$Symbol)
colnames(ieg15_hm) <- c("Control", "15_minute")
ieg60_hm <- cbind(ieg60$z0, ieg60$z15)
rownames(ieg60_hm) <- c(ieg60$Symbol)
colnames(ieg60_hm) <- c("Control", "60_minute")
ieg15UPDOWN_hm <- cbind(ieg15UPDOWN$z0, ieg15UPDOWN$z15)
rownames(ieg15UPDOWN_hm) <- c(ieg15UPDOWN$Symbol)
colnames(ieg15UPDOWN_hm) <- c("Control", "15_minute")
ieg60UPDOWN_hm <- cbind(ieg60UPDOWN$z0, ieg60UPDOWN$z15)
rownames(ieg60UPDOWN_hm) <- c(ieg60UPDOWN$Symbol)
colnames(ieg60UPDOWN_hm) <- c("Control", "60_minute")
################################################################################ piwi ##################################################3
piwi15_hm <- cbind(piwi15$z0, piwi15$z15)
rownames(piwi15_hm) <- c(piwi15$Transcript...Blast.Sequence.Features...Blast.Domain...Symbol)
colnames(piwi15_hm) <- c("Control", "15_minute")
piwi60_hm <- cbind(piwi60$z0, piwi60$z15)
rownames(piwi60_hm) <- c(piwi60$Transcript...Blast.Sequence.Features...Blast.Domain...Symbol)
colnames(piwi60_hm) <- c("Control", "60_minute")
piwi15UPDOWN_hm <- cbind(piwi15UPDOWN$z0, piwi15UPDOWN$z15)
rownames(piwi15UPDOWN_hm) <- c(piwi15UPDOWN$Transcript...Blast.Sequence.Features...Blast.Domain...Symbol)
colnames(piwi15UPDOWN_hm) <- c("Control", "15_minute")
piwi60UPDOWN_hm <- cbind(piwi60UPDOWN$z0, piwi60UPDOWN$z15)
rownames(piwi60UPDOWN_hm) <- c(piwi60UPDOWN$Transcript...Blast.Sequence.Features...Blast.Domain...Symbol)
colnames(piwi60UPDOWN_hm) <- c("Control", "60_minute")
################################################################################ NHEJ ##################################################3
NHEJ15_hm <- cbind(NHEJ15$z0, NHEJ15$z15)
rownames(NHEJ15_hm) <- c(NHEJ15$Transcript...Blast.Sequence.Features...Blast.Domain...Symbol)
colnames(NHEJ15_hm) <- c("Control", "15_minute")
NHEJ60_hm <- cbind(NHEJ60$z0, NHEJ60$z15)
rownames(NHEJ60_hm) <- c(NHEJ60$Transcript...Blast.Sequence.Features...Blast.Domain...Symbol)
colnames(NHEJ60_hm) <- c("Control", "60_minute")
NHEJ15UPDOWN_hm <- cbind(NHEJ15UPDOWN$z0, NHEJ15UPDOWN$z15)
rownames(NHEJ15UPDOWN_hm) <- c(NHEJ15UPDOWN$Transcript...Blast.Sequence.Features...Blast.Domain...Symbol)
colnames(NHEJ15UPDOWN_hm) <- c("Control", "15_minute")
NHEJ60UPDOWN_hm <- cbind(NHEJ60UPDOWN$z0, NHEJ60UPDOWN$z15)
rownames(NHEJ60UPDOWN_hm) <- c(NHEJ60UPDOWN$Transcript...Blast.Sequence.Features...Blast.Domain...Symbol)
colnames(NHEJ60UPDOWN_hm) <- c("Control", "60_minute")
################################################################################ NBFIG ##################################################3
NBFIG15_hm <- cbind(NBFIG15$z0, NBFIG15$z15)
rownames(NBFIG15_hm) <- c(NBFIG15$Marker)
colnames(NBFIG15_hm) <- c("Control", "15_minute")
NBFIG60_hm <- cbind(NBFIG60$z0, NBFIG60$z15)
rownames(NBFIG60_hm) <- c(NBFIG60$Marker)
colnames(NBFIG60_hm) <- c("Control", "60_minute")
NBFIG15UPDOWN_hm <- cbind(NBFIG15UPDOWN$z0, NBFIG15UPDOWN$z15)
rownames(NBFIG15UPDOWN_hm) <- c(NBFIG15UPDOWN$Marker)
colnames(NBFIG15UPDOWN_hm) <- c("Control", "15_minute")
NBFIG60UPDOWN_hm <- cbind(NBFIG60UPDOWN$z0, NBFIG60UPDOWN$z15)
rownames(NBFIG60UPDOWN_hm) <- c(NBFIG60UPDOWN$Marker)
colnames(NBFIG60UPDOWN_hm) <- c("Control", "60_minute")9 Plotting heatmaps of pathways showing significantly differentially expressed transcipts
##### Heatmaps of sigDE genes for each comparison plotting the
##### mean zscores####
Heatmap(calcium15_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(calcium15UPDOWN_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(calcium60_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(calcium60UPDOWN_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(cellcycle15_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(cellcycle15UPDOWN_hm, name = "Z-score", col = mycolz,
column_names_gp = gpar(fontsize = 8), cluster_columns = FALSE,
cluster_rows = TRUE)
Heatmap(cellcycle60_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(cellcycle60UPDOWN_hm, name = "Z-score", col = mycolz,
column_names_gp = gpar(fontsize = 8), cluster_columns = FALSE,
cluster_rows = TRUE)
Heatmap(celldeath15_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(celldeath15UPDOWN_hm, name = "Z-score", col = mycolz,
column_names_gp = gpar(fontsize = 8), cluster_columns = FALSE,
cluster_rows = TRUE)
Heatmap(celldeath60_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(celldeath60UPDOWN_hm, name = "Z-score", col = mycolz,
column_names_gp = gpar(fontsize = 8), cluster_columns = FALSE,
cluster_rows = TRUE)
Heatmap(cellmigration15_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(cellmigration15UPDOWN_hm, name = "Z-score", col = mycolz,
column_names_gp = gpar(fontsize = 8), cluster_columns = FALSE,
cluster_rows = TRUE)
Heatmap(cellmigration60_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(cellmigration60UPDOWN_hm, name = "Z-score", col = mycolz,
column_names_gp = gpar(fontsize = 8), cluster_columns = FALSE,
cluster_rows = TRUE)
Heatmap(DNAdamage15_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(DNAdamage15UPDOWN_hm, name = "Z-score", col = mycolz,
column_names_gp = gpar(fontsize = 8), cluster_columns = FALSE,
cluster_rows = TRUE)
Heatmap(DNAdamage60_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(DNAdamage60UPDOWN_hm, name = "Z-score", col = mycolz,
column_names_gp = gpar(fontsize = 8), cluster_columns = FALSE,
cluster_rows = TRUE)
Heatmap(DNArepair15_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(DNArepair15UPDOWN_hm, name = "Z-score", col = mycolz,
column_names_gp = gpar(fontsize = 8), cluster_columns = FALSE,
cluster_rows = TRUE)
Heatmap(DNArepair60_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(DNArepair60UPDOWN_hm, name = "Z-score", col = mycolz,
column_names_gp = gpar(fontsize = 8), cluster_columns = FALSE,
cluster_rows = TRUE)
Heatmap(HR15_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(HR15UPDOWN_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(HR60_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(HR60UPDOWN_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(ieg15_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(ieg15UPDOWN_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(ieg60_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(ieg60UPDOWN_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(MMR15_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(MMR15UPDOWN_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(MMR60_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(MMR60UPDOWN_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(NBFIG15_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE, split = NBFIG15$Cell.type)
Heatmap(NBFIG15UPDOWN_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE, split = NBFIG15UPDOWN$Cell.type)
Heatmap(NBFIG60_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE, split = NBFIG60$Cell.type)
Heatmap(NBFIG60UPDOWN_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE, split = NBFIG60UPDOWN$Cell.type)
Heatmap(nbsc15_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE, split = nbsc15$NB.subclass)
Heatmap(nbsc15UPDOWN_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE, split = nbsc15UPDOWN$NB.subclass)
Heatmap(nbsc60_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE, split = nbsc60$NB.subclass)
Heatmap(nbsc60UPDOWN_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE, split = nbsc60UPDOWN$NB.subclass)
Heatmap(neural15_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(neural15UPDOWN_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(neural60_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(neural60UPDOWN_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(NHEJ15_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(NHEJ15UPDOWN_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(NHEJ60_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(NHEJ60UPDOWN_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(piwi15_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(piwi15UPDOWN_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(piwi60_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(piwi60UPDOWN_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(proliferation15_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(proliferation15UPDOWN_hm, name = "Z-score", col = mycolz,
column_names_gp = gpar(fontsize = 8), cluster_columns = FALSE,
cluster_rows = TRUE)
Heatmap(proliferation60_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(proliferation60UPDOWN_hm, name = "Z-score", col = mycolz,
column_names_gp = gpar(fontsize = 8), cluster_columns = FALSE,
cluster_rows = TRUE)
Heatmap(replication15_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(replication15UPDOWN_hm, name = "Z-score", col = mycolz,
column_names_gp = gpar(fontsize = 8), cluster_columns = FALSE,
cluster_rows = TRUE)
Heatmap(replication60_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(replication60UPDOWN_hm, name = "Z-score", col = mycolz,
column_names_gp = gpar(fontsize = 8), cluster_columns = FALSE,
cluster_rows = TRUE)
Heatmap(signaling15_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(signaling15UPDOWN_hm, name = "Z-score", col = mycolz,
column_names_gp = gpar(fontsize = 8), cluster_columns = FALSE,
cluster_rows = TRUE)
Heatmap(signaling60_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE)
Heatmap(signaling60UPDOWN_hm, name = "Z-score", col = mycolz,
column_names_gp = gpar(fontsize = 8), cluster_columns = FALSE,
cluster_rows = TRUE)
Heatmap(sublethal15_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE, split = sublethal15$Sub.lethal..SL...cell.cluster)
Heatmap(sublethal15UPDOWN_hm, name = "Z-score", col = mycolz,
column_names_gp = gpar(fontsize = 8), cluster_columns = FALSE,
cluster_rows = TRUE, split = sublethal15UPDOWN$Sub.lethal..SL...cell.cluster)
Heatmap(sublethal60_hm, name = "Z-score", col = mycolz, column_names_gp = gpar(fontsize = 8),
cluster_columns = FALSE, cluster_rows = TRUE, split = sublethal60$Sub.lethal..SL...cell.cluster)
Heatmap(sublethal60UPDOWN_hm, name = "Z-score", col = mycolz,
column_names_gp = gpar(fontsize = 8), cluster_columns = FALSE,
cluster_rows = TRUE, split = sublethal60UPDOWN$Sub.lethal..SL...cell.cluster)
10 Gene set enrichment analysis for 15 minute and 60 minute timepoints
Next a gene enrichment analysis was performed using the topGO platform downloaded from Bioconductor here: :https://bioconductor.org/packages/release/bioc/html/topGO.html. The documentation for this platform can be found here: https://bioconductor.org/packages/release/bioc/vignettes/topGO/inst/doc/topGO.pdf. I used gene ontology terms annotated to the dd_Smed_v6 transcriptome from the planmine database found here: http://planmine.mpi-cbg.de/planmine/begin.do. Additionally I used the gene rankings from my analysis above to determine pathway enrichment. The Kolmogorv-Smirnov test is used to determine pathway significance.
Loading required package: BiocGenericsLoading required package: parallel
Attaching package: 'BiocGenerics'The following objects are masked from 'package:parallel':
clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
clusterExport, clusterMap, parApply, parCapply, parLapply,
parLapplyLB, parRapply, parSapply, parSapplyLBThe following objects are masked from 'package:dplyr':
combine, intersect, setdiff, unionThe following object is masked from 'package:limma':
plotMAThe following objects are masked from 'package:stats':
IQR, mad, sd, var, xtabsThe following objects are masked from 'package:base':
anyDuplicated, append, as.data.frame, basename, cbind, colnames,
dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
union, unique, unsplit, which.max, which.minLoading required package: graph
Attaching package: 'graph'The following object is masked from 'package:circlize':
degreeLoading required package: BiobaseWelcome to Bioconductor
Vignettes contain introductory material; view with
'browseVignettes()'. To cite Bioconductor, see
'citation("Biobase")', and for packages 'citation("pkgname")'.Loading required package: GO.dbLoading required package: AnnotationDbiLoading required package: stats4Loading required package: IRangesLoading required package: S4Vectors
Attaching package: 'S4Vectors'The following objects are masked from 'package:dplyr':
first, renameThe following object is masked from 'package:gplots':
spaceThe following object is masked from 'package:base':
expand.grid
Attaching package: 'IRanges'The following objects are masked from 'package:dplyr':
collapse, desc, slice
Attaching package: 'AnnotationDbi'The following object is masked from 'package:dplyr':
selectLoading required package: SparseM
Attaching package: 'SparseM'The following object is masked from 'package:base':
backsolve
groupGOTerms: GOBPTerm, GOMFTerm, GOCCTerm environments built.
Attaching package: 'topGO'The following object is masked from 'package:IRanges':
membersThe following object is masked from 'package:grid':
depth
Attaching package: 'genefilter'The following object is masked from 'package:ComplexHeatmap':
dist2genes <- read.delim("2018.5.7_geneID2GO.txt")
annot <- genes[, 1:2]
names(annot) <- c("gene", "GO_ID")
write.table(annot, file = "2018.5.20_gene2GO.map", sep = "\t",
quote = F, row.names = F, col.names = F)
geneID2GO <- readMappings(file = "2018.5.20_gene2GO.map")
head(geneID2GO)$dd_Smed_v6_10003_0_1
[1] "GO:0045454"
$dd_Smed_v6_100044_0_1
[1] "GO:0003777"
$dd_Smed_v6_100044_0_1
[1] "GO:0005524"
$dd_Smed_v6_100044_0_1
[1] "GO:0007018"
$dd_Smed_v6_100044_0_1
[1] "GO:0008017"
$dd_Smed_v6_10006_0_1
[1] "GO:0001522"ctrlv15 <- as.matrix(toptable_15vctrlred)
ctrlv60 <- as.matrix(toptable_60vctrlred)
ctrlv15FDR <- ctrlv15[, 5]
ctrlv60FDR <- ctrlv60[, 5]
topDiffGenes <- function(allScore) {
+return(allScore < 0.05)
}
##### 15 min v ctrl ############################# Create topGO
##### data object for Biological Process
GOdata <- new("topGOdata", description = "Gene Enrichment with 15 mins pDCS",
ontology = "BP", allGenes = ctrlv15FDR, geneSelectionFun = topDiffGenes,
annot = annFUN.gene2GO, gene2GO = geneID2GO)
Building most specific GOs ..... ( 176 GO terms found. )
Build GO DAG topology .......... ( 736 GO terms and 1472 relations. )
Annotating nodes ............... ( 871 genes annotated to the GO terms. )[1] 871[1] "dd_Smed_v6_663_0_1" "dd_Smed_v6_1595_0_1" "dd_Smed_v6_1769_0_1"
[4] "dd_Smed_v6_219_0_1" "dd_Smed_v6_8024_0_1" "dd_Smed_v6_257_0_1" A graphNEL graph with directed edges
Number of Nodes = 736
Number of Edges = 1472 # Creates list of significant gene IDs but not with
# associated GO terms
sg <- sigGenes(GOdata)
sg <- sigGenes(GOdata)
numSigGenes(GOdata)[1] 117write.table(sg, file = "2021.2.23_15vctrl_BP_GOgenelist.txt",
sep = "\t", quote = F, row.names = F, col.names = F)
# Kolmogorov-Smirnov testing
resultKS <- runTest(GOdata, algorithm = "elim", statistic = "ks")
-- Elim Algorithm --
the algorithm is scoring 736 nontrivial nodes
parameters:
test statistic: ks
cutOff: 0.01
score order: increasing
Level 14: 3 nodes to be scored (0 eliminated genes)
Level 13: 7 nodes to be scored (0 eliminated genes)
Level 12: 16 nodes to be scored (0 eliminated genes)
Level 11: 30 nodes to be scored (0 eliminated genes)
Level 10: 44 nodes to be scored (0 eliminated genes)
Level 9: 81 nodes to be scored (0 eliminated genes)
Level 8: 98 nodes to be scored (0 eliminated genes)
Level 7: 102 nodes to be scored (4 eliminated genes)
Level 6: 115 nodes to be scored (4 eliminated genes)
Level 5: 120 nodes to be scored (4 eliminated genes)
Level 4: 63 nodes to be scored (4 eliminated genes)
Level 3: 43 nodes to be scored (4 eliminated genes)
Level 2: 13 nodes to be scored (4 eliminated genes)
Level 1: 1 nodes to be scored (503 eliminated genes)tab <- GenTable(GOdata, KS = resultKS, topNodes = length(resultKS@score))
write.table(tab, file = "2021.2.23_15vctrl_BP_topnodes.txt",
sep = "\t", quote = F, row.names = F, col.names = F)
################################################# Create topGO data object for Cellular Component
GOdata <- new("topGOdata", description = "Gene Enrichment with 15 mins pDCS",
ontology = "CC", allGenes = ctrlv15FDR, geneSelectionFun = topDiffGenes,
annot = annFUN.gene2GO, gene2GO = geneID2GO)
Building most specific GOs ..... ( 94 GO terms found. )
Build GO DAG topology .......... ( 248 GO terms and 424 relations. )
Annotating nodes ............... ( 944 genes annotated to the GO terms. )[1] 944[1] "dd_Smed_v6_4005_0_1" "dd_Smed_v6_2696_0_1" "dd_Smed_v6_1089_0_1"
[4] "dd_Smed_v6_7837_0_1" "dd_Smed_v6_1731_0_1" "dd_Smed_v6_2119_0_1"A graphNEL graph with directed edges
Number of Nodes = 248
Number of Edges = 424 # Creates list of significant gene IDs but not with
# associated GO terms
sg <- sigGenes(GOdata)
sg <- sigGenes(GOdata)
numSigGenes(GOdata)[1] 125write.table(sg, file = "2021.2.23_15vctrl_CC_GOgenelist.txt",
sep = "\t", quote = F, row.names = F, col.names = F)
# Kolmogorov-Smirnov testing
resultKS <- runTest(GOdata, algorithm = "elim", statistic = "ks")
-- Elim Algorithm --
the algorithm is scoring 248 nontrivial nodes
parameters:
test statistic: ks
cutOff: 0.01
score order: increasing
Level 13: 1 nodes to be scored (0 eliminated genes)
Level 12: 5 nodes to be scored (0 eliminated genes)
Level 11: 8 nodes to be scored (0 eliminated genes)
Level 10: 16 nodes to be scored (0 eliminated genes)
Level 9: 21 nodes to be scored (0 eliminated genes)
Level 8: 34 nodes to be scored (0 eliminated genes)
Level 7: 41 nodes to be scored (5 eliminated genes)
Level 6: 32 nodes to be scored (5 eliminated genes)
Level 5: 30 nodes to be scored (5 eliminated genes)
Level 4: 32 nodes to be scored (5 eliminated genes)
Level 3: 23 nodes to be scored (16 eliminated genes)
Level 2: 4 nodes to be scored (16 eliminated genes)
Level 1: 1 nodes to be scored (16 eliminated genes)tab <- GenTable(GOdata, KS = resultKS, topNodes = length(resultKS@score))
write.table(tab, file = "2021.2.23_15vctrl_CC_topnodes.txt",
sep = "\t", quote = F, row.names = F, col.names = F)
################################################# Create topGO data object for Molecular Function
GOdata <- new("topGOdata", description = "Gene Enrichment with 15 mins pDCS",
ontology = "MF", allGenes = ctrlv15FDR, geneSelectionFun = topDiffGenes,
annot = annFUN.gene2GO, gene2GO = geneID2GO)
Building most specific GOs ..... ( 363 GO terms found. )
Build GO DAG topology .......... ( 692 GO terms and 903 relations. )
Annotating nodes ............... ( 7047 genes annotated to the GO terms. )[1] 7047[1] "dd_Smed_v6_659_0_1" "dd_Smed_v6_3281_0_1" "dd_Smed_v6_11100_0_1"
[4] "dd_Smed_v6_5635_0_1" "dd_Smed_v6_12472_0_1" "dd_Smed_v6_2075_0_1" A graphNEL graph with directed edges
Number of Nodes = 692
Number of Edges = 903 # Creates list of significant gene IDs but not with
# associated GO terms
sg <- sigGenes(GOdata)
sg <- sigGenes(GOdata)
numSigGenes(GOdata)[1] 914write.table(sg, file = "2021.2.23_15vctrl_MF_GOgenelist.txt",
sep = "\t", quote = F, row.names = F, col.names = F)
# Kolmogorov-Smirnov testing
resultKS <- runTest(GOdata, algorithm = "elim", statistic = "ks")
-- Elim Algorithm --
the algorithm is scoring 692 nontrivial nodes
parameters:
test statistic: ks
cutOff: 0.01
score order: increasing
Level 12: 1 nodes to be scored (0 eliminated genes)
Level 11: 3 nodes to be scored (0 eliminated genes)
Level 10: 14 nodes to be scored (0 eliminated genes)
Level 9: 34 nodes to be scored (0 eliminated genes)
Level 8: 62 nodes to be scored (0 eliminated genes)
Level 7: 103 nodes to be scored (0 eliminated genes)
Level 6: 190 nodes to be scored (0 eliminated genes)
Level 5: 142 nodes to be scored (36 eliminated genes)
Level 4: 94 nodes to be scored (167 eliminated genes)
Level 3: 37 nodes to be scored (167 eliminated genes)
Level 2: 11 nodes to be scored (1504 eliminated genes)
Level 1: 1 nodes to be scored (2041 eliminated genes)tab <- GenTable(GOdata, KS = resultKS, topNodes = length(resultKS@score))
write.table(tab, file = "2021.2.23_15vctrl_MF_topnodes.txt",
sep = "\t", quote = F, row.names = F, col.names = F)
##### 60 min v ctrl ############################# Create topGO
##### data object for Biological Process
GOdata <- new("topGOdata", description = "Gene Enrichment with 60 mins pDCS",
ontology = "BP", allGenes = ctrlv60FDR, geneSelectionFun = topDiffGenes,
annot = annFUN.gene2GO, gene2GO = geneID2GO)
Building most specific GOs ..... ( 176 GO terms found. )
Build GO DAG topology .......... ( 736 GO terms and 1472 relations. )
Annotating nodes ............... ( 871 genes annotated to the GO terms. )[1] 871[1] "dd_Smed_v6_3230_0_1" "dd_Smed_v6_1231_0_1" "dd_Smed_v6_1595_0_1"
[4] "dd_Smed_v6_147_0_1" "dd_Smed_v6_8163_0_1" "dd_Smed_v6_219_0_1" A graphNEL graph with directed edges
Number of Nodes = 736
Number of Edges = 1472 # Creates list of significant gene IDs but not with
# associated GO terms
sg <- sigGenes(GOdata)
sg <- sigGenes(GOdata)
numSigGenes(GOdata)[1] 121write.table(sg, file = "2021.2.23_60vctrl_BP_GOgenelist.txt",
sep = "\t", quote = F, row.names = F, col.names = F)
# Kolmogorov-Smirnov testing
resultKS <- runTest(GOdata, algorithm = "elim", statistic = "ks")
-- Elim Algorithm --
the algorithm is scoring 736 nontrivial nodes
parameters:
test statistic: ks
cutOff: 0.01
score order: increasing
Level 14: 3 nodes to be scored (0 eliminated genes)
Level 13: 7 nodes to be scored (0 eliminated genes)
Level 12: 16 nodes to be scored (0 eliminated genes)
Level 11: 30 nodes to be scored (0 eliminated genes)
Level 10: 44 nodes to be scored (0 eliminated genes)
Level 9: 81 nodes to be scored (0 eliminated genes)
Level 8: 98 nodes to be scored (0 eliminated genes)
Level 7: 102 nodes to be scored (4 eliminated genes)
Level 6: 115 nodes to be scored (4 eliminated genes)
Level 5: 120 nodes to be scored (4 eliminated genes)
Level 4: 63 nodes to be scored (4 eliminated genes)
Level 3: 43 nodes to be scored (8 eliminated genes)
Level 2: 13 nodes to be scored (9 eliminated genes)
Level 1: 1 nodes to be scored (9 eliminated genes)tab <- GenTable(GOdata, KS = resultKS, topNodes = length(resultKS@score))
write.table(tab, file = "2021.2.23_60vctrl_BP_topnodes.txt",
sep = "\t", quote = F, row.names = F, col.names = F)
################################################# Create topGO data object for Cellular Component
GOdata <- new("topGOdata", description = "Gene Enrichment with 60 mins pDCS",
ontology = "CC", allGenes = ctrlv60FDR, geneSelectionFun = topDiffGenes,
annot = annFUN.gene2GO, gene2GO = geneID2GO)
Building most specific GOs ..... ( 94 GO terms found. )
Build GO DAG topology .......... ( 248 GO terms and 424 relations. )
Annotating nodes ............... ( 944 genes annotated to the GO terms. )[1] 944[1] "dd_Smed_v6_2696_0_1" "dd_Smed_v6_1089_0_1" "dd_Smed_v6_4005_0_1"
[4] "dd_Smed_v6_2119_0_1" "dd_Smed_v6_6150_0_1" "dd_Smed_v6_2051_0_1"A graphNEL graph with directed edges
Number of Nodes = 248
Number of Edges = 424 # Creates list of significant gene IDs but not with
# associated GO terms
sg <- sigGenes(GOdata)
sg <- sigGenes(GOdata)
numSigGenes(GOdata)[1] 119write.table(sg, file = "2021.2.23_60vctrl_CC_GOgenelist.txt",
sep = "\t", quote = F, row.names = F, col.names = F)
# Kolmogorov-Smirnov testing
resultKS <- runTest(GOdata, algorithm = "elim", statistic = "ks")
-- Elim Algorithm --
the algorithm is scoring 248 nontrivial nodes
parameters:
test statistic: ks
cutOff: 0.01
score order: increasing
Level 13: 1 nodes to be scored (0 eliminated genes)
Level 12: 5 nodes to be scored (0 eliminated genes)
Level 11: 8 nodes to be scored (0 eliminated genes)
Level 10: 16 nodes to be scored (0 eliminated genes)
Level 9: 21 nodes to be scored (0 eliminated genes)
Level 8: 34 nodes to be scored (0 eliminated genes)
Level 7: 41 nodes to be scored (0 eliminated genes)
Level 6: 32 nodes to be scored (0 eliminated genes)
Level 5: 30 nodes to be scored (103 eliminated genes)
Level 4: 32 nodes to be scored (103 eliminated genes)
Level 3: 23 nodes to be scored (103 eliminated genes)
Level 2: 4 nodes to be scored (103 eliminated genes)
Level 1: 1 nodes to be scored (103 eliminated genes)tab <- GenTable(GOdata, KS = resultKS, topNodes = length(resultKS@score))
write.table(tab, file = "2021.2.23_60vctrl_CC_topnodes.txt",
sep = "\t", quote = F, row.names = F, col.names = F)
################################################# Create topGO data object for Molecular Function
GOdata <- new("topGOdata", description = "Gene Enrichment with 60 mins pDCS",
ontology = "MF", allGenes = ctrlv60FDR, geneSelectionFun = topDiffGenes,
annot = annFUN.gene2GO, gene2GO = geneID2GO)
Building most specific GOs ..... ( 363 GO terms found. )
Build GO DAG topology .......... ( 692 GO terms and 903 relations. )
Annotating nodes ............... ( 7047 genes annotated to the GO terms. )[1] 7047[1] "dd_Smed_v6_4392_0_1" "dd_Smed_v6_7295_0_1" "dd_Smed_v6_7144_0_1"
[4] "dd_Smed_v6_11100_0_1" "dd_Smed_v6_3281_0_1" "dd_Smed_v6_1913_0_1" A graphNEL graph with directed edges
Number of Nodes = 692
Number of Edges = 903 # Creates list of significant gene IDs but not with
# associated GO terms
sg <- sigGenes(GOdata)
sg <- sigGenes(GOdata)
numSigGenes(GOdata)[1] 893write.table(sg, file = "2021.2.23_60vctrl_MF_GOgenelist.txt",
sep = "\t", quote = F, row.names = F, col.names = F)
# Kolmogorov-Smirnov testing
resultKS <- runTest(GOdata, algorithm = "elim", statistic = "ks")
-- Elim Algorithm --
the algorithm is scoring 692 nontrivial nodes
parameters:
test statistic: ks
cutOff: 0.01
score order: increasing
Level 12: 1 nodes to be scored (0 eliminated genes)
Level 11: 3 nodes to be scored (0 eliminated genes)
Level 10: 14 nodes to be scored (0 eliminated genes)
Level 9: 34 nodes to be scored (0 eliminated genes)
Level 8: 62 nodes to be scored (0 eliminated genes)
Level 7: 103 nodes to be scored (0 eliminated genes)
Level 6: 190 nodes to be scored (25 eliminated genes)
Level 5: 142 nodes to be scored (296 eliminated genes)
Level 4: 94 nodes to be scored (751 eliminated genes)
Level 3: 37 nodes to be scored (751 eliminated genes)
Level 2: 11 nodes to be scored (874 eliminated genes)
Level 1: 1 nodes to be scored (2645 eliminated genes)tab <- GenTable(GOdata, KS = resultKS, topNodes = length(resultKS@score))
write.table(tab, file = "2021.2.23_60vctrl_MF_topnodes.txt",
sep = "\t", quote = F, row.names = F, col.names = F)
# looking at termination of G-protein coupled recepto... GO
# term goID <- tab2[1, 'GO.ID']
# print(showGroupDensity(GOdata, goID, ranks = TRUE))#Session Info
R version 4.0.3 (2020-10-10)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur 10.16
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats4 parallel grid stats graphics grDevices utils
[8] datasets methods base
other attached packages:
[1] genefilter_1.72.1 topGO_2.42.0 SparseM_1.81
[4] GO.db_3.12.1 AnnotationDbi_1.52.0 IRanges_2.24.1
[7] S4Vectors_0.28.1 Biobase_2.50.0 graph_1.68.0
[10] BiocGenerics_0.36.0 circlize_0.4.12 ComplexHeatmap_2.6.2
[13] RColorBrewer_1.1-2 dplyr_1.0.4 edgeR_3.32.1
[16] limma_3.46.0 ggplot2_3.3.3 gplots_3.1.1
loaded via a namespace (and not attached):
[1] httr_1.4.2 sass_0.3.1 splines_4.0.3
[4] bit64_4.0.5 jsonlite_1.7.2 gtools_3.8.2
[7] bslib_0.2.4 assertthat_0.2.1 highr_0.8
[10] blob_1.2.1 yaml_2.2.1 pillar_1.5.0
[13] RSQLite_2.2.3 lattice_0.20-41 glue_1.4.2
[16] digest_0.6.27 colorspace_2.0-0 Matrix_1.3-2
[19] htmltools_0.5.1.1 XML_3.99-0.5 pkgconfig_2.0.3
[22] GetoptLong_1.0.5 xtable_1.8-4 purrr_0.3.4
[25] scales_1.1.1 annotate_1.68.0 tibble_3.0.6
[28] generics_0.1.0 ellipsis_0.3.1 cachem_1.0.4
[31] withr_2.4.1 survival_3.2-7 magrittr_2.0.1
[34] crayon_1.4.1 memoise_2.0.0 evaluate_0.14
[37] fansi_0.4.2 Cairo_1.5-12.2 tools_4.0.3
[40] GlobalOptions_0.1.2 formatR_1.7 lifecycle_1.0.0
[43] matrixStats_0.58.0 stringr_1.4.0 munsell_0.5.0
[46] locfit_1.5-9.4 cluster_2.1.1 compiler_4.0.3
[49] jquerylib_0.1.3 caTools_1.18.1 rlang_0.4.10
[52] rjson_0.2.20 bitops_1.0-6 rmarkdown_2.7
[55] gtable_0.3.0 DBI_1.1.1 R6_2.5.0
[58] knitr_1.31 fastmap_1.1.0 bit_4.0.4
[61] utf8_1.1.4 clue_0.3-58 KernSmooth_2.23-18
[64] shape_1.4.5 stringi_1.5.3 Rcpp_1.0.6
[67] vctrs_0.3.6 png_0.1-7 tidyselect_1.1.0
[70] xfun_0.22